با استفاده از داده های زلزله ها در ایران و جهان به سوالات زیر پاسخ دهید.
library(readr)
library(dplyr)
library(ggplot2)
library(plotly)
library(highcharter)
library(raster)
library(sp)
library(RColorBrewer)
library(ggmap)
library(ggthemes)
historical_web_data = read_rds("D:/University/Semester8/Data analysis/week11/data/historical_web_data_26112015.rds")
disaster = read_delim("D:/University/Semester8/Data analysis/week11/data/disaster.txt", "\t", escape_double = FALSE, trim_ws = TRUE)
worldwide = read.csv("D:/University/Semester8/Data analysis/week11/data/worldwide.csv")
iran_earthquake = read_rds("D:/University/Semester8/Data analysis/week11/data/iran_earthquake.rds")۱. با استفاده از داده های historical_web_data_26112015.rds و استفاده از نمودار پراکنش سه بعدی بسته plotly نمودار طول، عرض و عمق زلزله ها را رسم نمایید. علاوه بر آن بزرگی هر نقطه را برابر بزرگی زمین لرزه قرار دهید.
plot_ly(historical_web_data , x = ~historical_web_data$Longitude , y = ~historical_web_data$Latitude , z = ~historical_web_data$Depth , size = ~historical_web_data$Magnitude)۲. پویانمایی سونامی های تاریخی را بر حسب شدت بر روی نقشه زمین رسم نمایید.(از داده زلزله های بزرگ استفاده نمایید.)
disaster %>% rename(lat = LATITUDE,lon = LONGITUDE, z = DEATHS,name = COUNTRY,sequence = YEAR) %>% dplyr::select(lat, lon, z, name, sequence) -> dis
hcmap() %>% hc_add_series(data = dis, type = "mapbubble",minSize = 0, maxSize = 30) %>% hc_plotOptions(series = list(showInLegend = FALSE))۳. نمودار چگالی دو بعدی زلزله های تاریخی ایران را رسم کنید.( از داده iran_earthquake.rds و لایه stat_density_2d استفاده نمایید).
iran_earthquake %>% filter(Mag > 5) -> iran_earthquake1
Iran <- map_data("world") %>% filter(region=="Iran")
iran_earthquake1 = iran_earthquake1 %>% filter(Mag > 5)
ggplot(iran_earthquake1 , aes(x = iran_earthquake1$Long , y = iran_earthquake1$Lat)) + geom_polygon(data = Iran , aes(x = Iran$long , y = Iran$lat)) + geom_point(color="red") + geom_density_2d() + xlab("Longitude") + ylab("Latitude")
۴. احتمال اینکه در ایران در پنج سال آینده زلزله به بزرگی هفت ریشتر رخ دهد را محاسبه کنید. (از احتمال شرطی استفاده کنید.)
powerful_earthquake = disaster %>% filter(COUNTRY == "IRAN" ,EQ_PRIMARY > 7)
y = powerful_earthquake $ YEAR
y = y[-1]
y1 = c(y[-1] , 0)
y = (y1 - y)[-26]
print("Probability that we have an earthauake more than 7 rishter for next 5 years if there has not been any since last one(2017) is(we are in 2018):")## [1] "Probability that we have an earthauake more than 7 rishter for next 5 years if there has not been any since last one(2017) is(we are in 2018):"print(sum(y<= 6)/length(y))## [1] 0.28۵. بر اساس داده های زلزله های بزرگ ابتدا تعداد و متوسط کشته زلزله ها را بر حسب کشور استخراج نمایید. سپس نمودار گرمایی تعداد کشته ها را بر روی کره زمین رسم نمایید.(مانند مثال زیر!)
tlb = NULL
for(i in unique(disaster$COUNTRY)){
x = disaster %>% filter(COUNTRY == i)
d = sum(x$DEATHS , na.rm = T)
ave = mean(x$DEATHS , na.rm = T)
x = c(i , d , ave)
tlb = rbind(tlb , x)
}
colnames(tlb) = c("region" , "DEATH" , "DEATH_AVERAGE")
tlb = as.data.frame(tlb)
world <- map_data("world")
world = world %>% filter(region != "Antarctica")
world$region = world$region %>% toupper()
a = unique(world$region)
b = unique(tlb$region)
c = b %in% a
b = b[c]
tlb = tlb %>% filter(region %in% b)
e = world %>% filter(region %in% b)
e1 = merge(x = e , y = tlb , by = "region" , all = T)
rr = as.numeric(as.character(e1$DEATH))
rr = cut(rr, quantile(rr, seq(0, 1, len = 10)),include.lowest = TRUE)
e2 = e %>% mutate(DEATHS = rr)
ggplot(e2) + geom_polygon(aes(long, lat, group = group , fill = DEATHS)) + theme_map()
***
۶. با استفاده از داده لرزه های بزرگ و به وسیله طول، عرض، شدت، عمق مدلی برای پیش بینی تعداد کشته های زلزله بیابید.
fit = glm(data = disaster , DEATHS ~ LATITUDE + LONGITUDE + FOCAL_DEPTH + EQ_PRIMARY)۷. با استفاده از داده worldwide.csv به چند سوال زیر پاسخ دهید. تحقیق کنید آیا می توان از پیش لرزه، زلزله اصلی را پیش بینی کرد؟
۸. گزاره " آیا شدت زلزله به عمق آن بستگی دارد" را تحقیق کنید؟ (طبیعتا از آزمون فرض باید استفاده کنید.)
chisq.test(worldwide$mag , worldwide$depth)##
## Pearson's Chi-squared test
##
## data: worldwide$mag and worldwide$depth
## X-squared = 3403700, df = 4463700, p-value = 1با توجه به p_value بالایی که به دست آوردیم این فرض که عمق زلزله از بزرگی مستقل است را رد نمی کنیم.
۹. میانگین سالانه زلزله ها را بر حسب کشور به دست آورید. آیا میتوان دلیلی در تایید یا رد تئوری هارپ ارائه کرد.
worldwide_by_year = split(worldwide, format(as.Date(worldwide$time), "%Y"))
loc = NULL
k = worldwide$place
k = as.character(k)
p = NULL
for(i in k){
x = strsplit(i , "[,]")
x = x[[1]][2]
p = c(p ,x)
}
p = substring(p , 2)
p1 = unique(p)
USA = p1[p1 %in% state.name]
USA = c(USA , "CA" , "U.S. Virgin Islands")
for(i in p1){
if(i %in% USA){
loc = c(loc , "USA")
}
else{
loc = c(loc , i)
}
}
loc = unique(loc)
tlb_year = as.data.frame(loc)
colnames(tlb_year) = "country"
for(j in worldwide_by_year){
p = NULL
k = j$place
k = as.character(k)
for(i in k){
x = strsplit(i , "[,]")
x = x[[1]][2]
p = c(p ,x)
}
p = substring(p , 2)
p1 = unique(p)
USA = p1[p1 %in% state.name]
USA = c(USA , "CA" , "U.S. Virgin Islands")
country = NULL
for(i in p){
if(i %in% USA){
country = c(country , "USA")
}
else{
country = c(country , i)
}
}
a = as.data.frame(table(country))
tlb_year = merge(x = tlb_year , y = a , by = "country" , all = T)
}
colnames(tlb_year) = c("Country" , "2015" , "2016" , "2017" , "2018")
tlb_year[is.na(tlb_year)] <- 0
tlb_year = tlb_year %>% mutate(mean_by_year = rowMeans(subset(tlb_year , select = c(`2015`, `2016`,`2017`,`2018`))))
print(tlb_year)## Country 2015 2016 2017 2018
## 1 Afghanistan 15 193 159 46
## 2 Albania 0 7 8 3
## 3 Algeria 0 15 7 4
## 4 American Samoa 0 1 1 1
## 5 Angola 0 1 0 0
## 6 Anguilla 0 15 11 1
## 7 Antarctica 0 5 5 2
## 8 Antigua and Barbuda 0 1 2 1
## 9 Argentina 18 234 261 61
## 10 Armenia 0 1 1 0
## 11 Australia 1 33 21 2
## 12 Austria 0 1 2 2
## 13 Azerbaijan 1 9 10 1
## 14 B.C. 10 87 64 22
## 15 Bangladesh 1 2 2 1
## 16 Barbados 1 1 4 0
## 17 Barbuda 0 6 10 1
## 18 Bhutan 0 1 2 0
## 19 Bolivia 2 45 43 22
## 20 Bonaire 0 3 3 2
## 21 Bosnia and Herzegovina 0 3 8 0
## 22 Botswana 0 0 12 0
## 23 Bouvet Island 1 5 3 1
## 24 Brazil 3 1 2 2
## 25 British Indian Ocean Territory 0 1 1 0
## 26 British Virgin Islands 32 674 423 100
## 27 Bulgaria 0 2 4 2
## 28 Burma 1 41 39 47
## 29 Canada 8 127 130 58
## 30 Cape Verde 0 0 1 0
## 31 Cayman Islands 1 0 3 1
## 32 Chile 103 776 764 160
## 33 China 12 204 165 53
## 34 Christmas Island 3 4 7 2
## 35 Colombia 3 106 82 28
## 36 Comoros 0 1 0 0
## 37 Costa Rica 2 33 32 19
## 38 Croatia 0 4 2 1
## 39 Cuba 0 13 8 1
## 40 Curaçao 0 0 1 0
## 41 Cyprus 2 12 7 0
## 42 Democratic Republic of the Congo 1 4 7 7
## 43 Djibouti 0 2 4 0
## 44 Dominica 0 6 1 0
## 45 Dominican Republic 24 202 252 194
## 46 East Timor 7 79 60 20
## 47 Ecuador 3 176 59 23
## 48 Ecuador region 0 4 1 0
## 49 Egypt 0 3 2 1
## 50 El Salvador 4 107 69 43
## 51 Equatorial Guinea 0 1 0 0
## 52 Eritrea 0 1 0 0
## 53 Ethiopia 0 1 1 1
## 54 Falkland Islands 0 1 2 0
## 55 Federated States of Micronesia 0 1 1 0
## 56 Fiji 40 630 621 151
## 57 Finland 0 0 1 0
## 58 France 1 15 8 6
## 59 French Guiana 0 0 2 0
## 60 French Polynesia 0 0 1 0
## 61 French Polynesia region 0 2 0 0
## 62 French Southern Territories 0 0 2 0
## 63 Germany 0 0 0 1
## 64 Greece 17 164 171 65
## 65 Greenland 1 8 10 11
## 66 Grenada 0 0 2 0
## 67 Guadeloupe 1 9 3 3
## 68 Guam 2 119 102 29
## 69 Guatemala 12 94 92 32
## 70 Guinea 0 0 0 1
## 71 Haiti 1 7 4 3
## 72 Honduras 1 17 7 22
## 73 Iceland 0 3 23 9
## 74 India 10 94 100 31
## 75 India region 0 0 2 1
## 76 Indonesia 244 1659 1523 464
## 77 Iran 6 110 204 87
## 78 Iraq 0 4 17 38
## 79 Italy 5 93 34 13
## 80 Jamaica 0 3 1 1
## 81 Japan 66 1165 805 239
## 82 Japan region 3 43 30 17
## 83 Kazakhstan 1 6 10 7
## 84 Kiribati region 0 1 0 1
## 85 Kosovo 0 0 0 1
## 86 Kyrgyzstan 3 83 27 8
## 87 Laikit II (Dimembe) 0 1 1 1
## 88 Laos 0 1 0 0
## 89 Lebanon 0 1 0 0
## 90 Lesotho 0 2 0 0
## 91 Libya 0 0 2 0
## 92 Macedonia 0 10 14 0
## 93 Madagascar 0 0 1 1
## 94 Malawi 0 2 3 6
## 95 Malaysia 0 1 1 2
## 96 Maldives 0 1 0 0
## 97 Malta 0 1 1 0
## 98 Martinique 1 8 3 1
## 99 Mauritania 1 1 0 0
## 100 Mauritius 0 6 6 1
## 101 Mexico 18 279 779 131
## 102 Micronesia 3 58 49 6
## 103 Mongolia 1 5 4 5
## 104 Montenegro 0 3 0 2
## 105 Montserrat 0 0 1 0
## 106 Morocco 0 48 2 0
## 107 Mozambique 0 3 3 1
## 108 Namibia 0 2 0 1
## 109 Nepal 2 25 9 0
## 110 New Caledonia 2 145 415 25
## 111 New Zealand 48 1076 443 120
## 112 Nicaragua 6 99 59 29
## 113 Niger 0 0 1 0
## 114 Niue 0 0 3 0
## 115 Norfolk Island 0 1 0 0
## 116 North Korea 0 2 6 0
## 117 Northern Mariana Islands 15 440 262 84
## 118 Northwest Territories 0 1 0 0
## 119 Norway 0 1 2 0
## 120 Nunavut 0 0 2 0
## 121 Nunoa 0 1 3 1
## 122 Oman 0 1 3 2
## 123 Pakistan 1 40 32 11
## 124 Palau 0 2 1 2
## 125 Panama 2 30 34 13
## 126 Papua New Guinea 60 946 607 529
## 127 Peru 27 281 231 65
## 128 Philippines 47 577 498 139
## 129 Poland 0 6 5 1
## 130 Portugal 0 20 15 7
## 131 Prince Edward Islands 0 6 0 1
## 132 Puerto Rico 104 1043 836 346
## 133 Reunion 0 0 0 1
## 134 Romania 2 5 12 4
## 135 Russia 25 365 482 124
## 136 Rwanda 0 3 0 0
## 137 Saint Barthélemy 0 1 1 2
## 138 Saint Helena 0 9 6 3
## 139 Saint Kitts and Nevis 0 1 3 2
## 140 Saint Lucia 0 0 1 0
## 141 Saint Martin 1 5 0 0
## 142 Saint Vincent and the Grenadines 0 0 1 0
## 143 Samoa 1 7 3 0
## 144 Saudi Arabia 0 8 3 1
## 145 Serbia 0 1 0 0
## 146 Solomon Islands 14 505 286 61
## 147 Somalia 0 10 7 3
## 148 South Africa 0 12 8 10
## 149 South Georgia and the South Sandwich Islands 14 331 174 41
## 150 South Korea 1 5 2 1
## 151 South Sandwich Islands 10 65 81 24
## 152 Spain 0 4 5 2
## 153 Sudan 0 0 1 0
## 154 Svalbard and Jan Mayen 1 17 17 4
## 155 Sweden 0 1 0 0
## 156 Switzerland 0 2 2 0
## 157 Syria 0 2 0 0
## 158 Taiwan 12 105 56 69
## 159 Tajikistan 41 146 111 49
## 160 Tanzania 0 8 10 4
## 161 Thailand 0 0 3 1
## 162 Tonga 36 442 432 108
## 163 Trinidad and Tobago 0 8 8 1
## 164 Tunisia 0 1 1 1
## 165 Turkey 6 42 89 16
## 166 Turkmenistan 2 5 19 1
## 167 Uganda 0 1 4 0
## 168 Ukraine 0 1 0 0
## 169 United Kingdom 0 0 2 2
## 170 USA 673 7825 7895 5775
## 171 Uzbekistan 1 10 5 3
## 172 Vanuatu 32 553 271 73
## 173 Venezuela 2 19 19 13
## 174 Vietnam 0 1 1 2
## 175 Wallis and Futuna 10 31 41 16
## 176 Yemen 0 19 17 3
## 177 Zambia 0 2 7 1
## 178 Zimbabwe 0 4 2 1
## 179 <NA> 0 0 0 0
## mean_by_year
## 1 103.25
## 2 4.50
## 3 6.50
## 4 0.75
## 5 0.25
## 6 6.75
## 7 3.00
## 8 1.00
## 9 143.50
## 10 0.50
## 11 14.25
## 12 1.25
## 13 5.25
## 14 45.75
## 15 1.50
## 16 1.50
## 17 4.25
## 18 0.75
## 19 28.00
## 20 2.00
## 21 2.75
## 22 3.00
## 23 2.50
## 24 2.00
## 25 0.50
## 26 307.25
## 27 2.00
## 28 32.00
## 29 80.75
## 30 0.25
## 31 1.25
## 32 450.75
## 33 108.50
## 34 4.00
## 35 54.75
## 36 0.25
## 37 21.50
## 38 1.75
## 39 5.50
## 40 0.25
## 41 5.25
## 42 4.75
## 43 1.50
## 44 1.75
## 45 168.00
## 46 41.50
## 47 65.25
## 48 1.25
## 49 1.50
## 50 55.75
## 51 0.25
## 52 0.25
## 53 0.75
## 54 0.75
## 55 0.50
## 56 360.50
## 57 0.25
## 58 7.50
## 59 0.50
## 60 0.25
## 61 0.50
## 62 0.50
## 63 0.25
## 64 104.25
## 65 7.50
## 66 0.50
## 67 4.00
## 68 63.00
## 69 57.50
## 70 0.25
## 71 3.75
## 72 11.75
## 73 8.75
## 74 58.75
## 75 0.75
## 76 972.50
## 77 101.75
## 78 14.75
## 79 36.25
## 80 1.25
## 81 568.75
## 82 23.25
## 83 6.00
## 84 0.50
## 85 0.25
## 86 30.25
## 87 0.75
## 88 0.25
## 89 0.25
## 90 0.50
## 91 0.50
## 92 6.00
## 93 0.50
## 94 2.75
## 95 1.00
## 96 0.25
## 97 0.50
## 98 3.25
## 99 0.50
## 100 3.25
## 101 301.75
## 102 29.00
## 103 3.75
## 104 1.25
## 105 0.25
## 106 12.50
## 107 1.75
## 108 0.75
## 109 9.00
## 110 146.75
## 111 421.75
## 112 48.25
## 113 0.25
## 114 0.75
## 115 0.25
## 116 2.00
## 117 200.25
## 118 0.25
## 119 0.75
## 120 0.50
## 121 1.25
## 122 1.50
## 123 21.00
## 124 1.25
## 125 19.75
## 126 535.50
## 127 151.00
## 128 315.25
## 129 3.00
## 130 10.50
## 131 1.75
## 132 582.25
## 133 0.25
## 134 5.75
## 135 249.00
## 136 0.75
## 137 1.00
## 138 4.50
## 139 1.50
## 140 0.25
## 141 1.50
## 142 0.25
## 143 2.75
## 144 3.00
## 145 0.25
## 146 216.50
## 147 5.00
## 148 7.50
## 149 140.00
## 150 2.25
## 151 45.00
## 152 2.75
## 153 0.25
## 154 9.75
## 155 0.25
## 156 1.00
## 157 0.50
## 158 60.50
## 159 86.75
## 160 5.50
## 161 1.00
## 162 254.50
## 163 4.25
## 164 0.75
## 165 38.25
## 166 6.75
## 167 1.25
## 168 0.25
## 169 1.00
## 170 5542.00
## 171 4.75
## 172 232.25
## 173 13.25
## 174 1.00
## 175 24.50
## 176 9.75
## 177 2.50
## 178 1.75
## 179 0.00۱۰. سه حقیقت جالب در مورد زلزله بیابید.
تعداد کشته های زلزله های بالای 7 ریشتر را برای سال های یعد از 1975 کشیدیم که سیر نزولی داشت. می توان گفت با پیشرفت تکنولوژی تعداد کشته های زازله ها کم شده است.
deaths_per_year = NULL
y = unique(disaster$YEAR)
disaster1 = disaster %>% filter(EQ_PRIMARY > 7)
for(i in y){
x = disaster1 %>% filter(YEAR == i)
x = sum(x$DEATHS , na.rm = T)
deaths_per_year = c(deaths_per_year , x)
}
k = (deaths_per_year != 0)
deaths_per_year = deaths_per_year[k]
y = y[k]
k = y>1975
y = y[k]
deaths_per_year = deaths_per_year[k]
tlb_death_per_year = data.frame(y , deaths_per_year)
colnames(tlb_death_per_year) = c("year" , "deaths")
ggplot(tlb_death_per_year , aes(y , deaths_per_year)) + geom_col() + geom_smooth(method = "lm", se = T)
توزیع ریشتر زلزله ها را از داده disaster کشیدیم که به یک توزیع نرمال با میانگین 6.47 رسیدیم.
y = disaster$EQ_PRIMARY
y <- y[!is.na(y)]
m = mean(y)
y = as.data.frame(table(y))
colnames(y) = c("Year" , "Freq")
ggplot(y , aes(Year , Freq)) + geom_col() + coord_flip()
نمودار 30 کشوری که بیشترین تعداد زلزله هر را دارا بودند بر حسب تعداد زلزله ها کشیدیم. مشاهده می شود تعداد زلزله های این 30 کشور حدود 0.8 کل زلزله ها را شامل می شود و کشور چین به تنهایی 0.1 تعداد کل زلزله ها را دارا است.
y = disaster$COUNTRY
y <- y[!is.na(y)]
m = mean(y)
y = as.data.frame(table(y))
y <- y[order(-y$Freq),]
ratio = sum(y$Freq[1:30]) / sum(y$Freq)
print(paste("30 maximim earthquakes number to sum of earthquakes : " , ratio , sep = ""))## [1] "30 maximim earthquakes number to sum of earthquakes : 0.806173249253236"ratio_china = y$Freq[1] / sum(y$Freq)
print(paste("China earthquakes number to sum of earthquakes : " , ratio_china , sep = ""))## [1] "China earthquakes number to sum of earthquakes : 0.0970793229339529"y = y[1:30,]
colnames(y) = c("Year" , "Freq")
print(ggplot(y , aes(reorder(Year,-Freq) , Freq)) + geom_col() + coord_flip())